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A DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ 
EQUATION BASED ON RAY THEORY 

CHRISTIAAN C. STOLK 


Abstract. We develop a new dispersion minimizing compact finite difference scheme for 
the Helmholtz equation in 2 and 3 dimensions. The scheme is based on a newly developed 
ray theory for difference equations. A discrete Helmholtz operator and a discrete operator 
to be applied to the source and the wavefields are constructed. Their coefficients are 
piecewise polynomial functions of hk, chosen such that phase and amplitude errors are 
minimal. The phase errors of the scheme are very small, approximately as small as 
those of the 2-D quasi-stabilized FEM method and substantially smaller than those of 
alternatives in 3-D, assuming the same number of gridpoints per wavelength is used. In 
numerical experiments, accurate solutions are obtained in constant and smoothly varying 
media using meshes with only five to six points per wavelength and wave propagation 
over hundreds of wavelengths. When used as a coarse level discretization in a multigrid 
method the scheme can even be used with downto three points per wavelength. Tests 
on 3-D examples with up to 10 8 degrees of freedom show that with a recently developed 
hybrid solver, the use of coarser meshes can lead to corresponding savings in computation 
time, resulting in good simulation times compared to the literature. 


1. Introduction 

We consider the discretization on regular meshes of the Helmholtz equation 
(1) — A u — k{x) 2 u — f 

with large and variable k. These methods are widely used for simulations on unbounded 
domains, for example in exploration geophysics, using domain sizes, in three dimensions, 
of up to hundreds of wavelengths [23; 3 E2 EH]- 

A key issue for such discretizations are the dispersion (phase) errors, that are closely 
related to pollution errors |3j. Typically, the propagating wave solutions to the discrete 
and continuous equations have slightly different wavelengths. These wavelength errors are 
also referred to as phase velocity or phase slowness errors, in which case they are differently 
normalized. They lead to phase errors in the solution that grow with the distance from 
the source. A second important consideration is solver cost. The discretized Helmholtz 
operator should of course be cheap to apply and/or invert. 

A class of discretizations, that performs relatively well on these criteria, is given by so 
called compact finite difference methods, that use a 3 x 3 square or 3 x 3 x 3 cubic stencil 
in two resp. three dimensions. The corresponding discrete Helmholtz operators can be 
efficiently applied and inverted compared for example to standard finite difference or finite 
element methods. Many authors have studied such discretizations and obtained formulae 
for the coefficients as a function of k and the grid spacing h [3, [13, Iff. 21. 2?, [6. 29, [26] . 
We will discuss these schemes more in detail below, and compare their phase slowness 
errors with those of standard finite differences and Lagrange finite element methods on 
regular meshes. 

To design such methods, several strategies have been followed. One approach is too 
construct schemes of higher order, for example order four order six, see mmm and the 
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references in m- Another approach is to stay with second order schemes but minimize the 
dispersion errors, because these are the dominant errors for long distance wave propagation 
mmmmm- From the point of view of phase slowness errors, the sixth order schemes 
of |27j and [29] and the quasi-stabilized FEM (QS-FEM) scheme of [3] are the best, see the 
results below. The latter has the smallest phase slowness errors by a substantial margin, 
but is only available in 2-D. 

Alternative methods include higher order finite elements. An advantage of these meth¬ 
ods is the better theory for the behavior of the errors in the limit that the grid spacing 
goes to zero, see for example mm- 

In this paper we introduce a new second order dispersion minimizing scheme in 2 and 3 
dimensions with phase slowness errors comparable to those of QS-FEM. Accurate ampli¬ 
tudes are obtained as well using new amplitude correction operators. A theoretical justifi¬ 
cation is given using a newly developed ray theory for Helmholtz-like difference equations. 
This theory is remarkably similar to the continuous theory, when both are formulated in 
terms of the symbols associated with the operators. With numerical examples we show 
the potential for accurate and fast simulation on relatively coarse meshes. In addition we 
show applications where the method is used as a coarse level discretization in multigrid 
solvers 

We will briefly describe the methodology and the results. It is known that the second 
order, compact finite difference discretizations of the Helmholtz operator form a 3 or 5 
parameter family, in 2 and 3 dimensions respectively, and that by choosing parameters in 
a certain way, the phase slowness errors can be reduced compared to standard schemes m 
HQ- When coefficients are allowed to depend on hk in a piecewise constant [6] or piecewise 
linear fashion [26j, they can be further reduced. In this paper we let the parameters 
depend in a C 1 fashion on hk through third order Hermite interpolation and obtain a 
further reduction of the phase slowness errors. 

Dispersion minimizing schemes are typically intended for use on quite coarse meshes, 
and a theoretical understanding that does not involve the limit hk -A 0 is therefore of 
considerable interest. For this reason we consider ray theory for Helmholtz-like difference 
equations. 

Ray theory for continuous Helmholtz equations is well known m- Solutions are sought 
in the form A(x)e luJ ®( x \ If k is smooth, $ satisfies a certain eikonal equation and A a 
certain transport equation, than such solutions approximate the true solutions increasingly 
well in the limit uo -A oo. Here we develop a similar theory for Helmholtz-like difference 
equations. We can then choose the discrete scheme such that the phase and amplitude 
functions associated with the discrete operator approximate match those of the continuous 
operator well. As can be expected, schemes with small phase slowness errors have accurate 
phase functions. By introducing amplitude correction operators, accurate amplitude of the 
ray-theoretic solutions are obtained. 

We are interested in two ways of applying the discretized Helmholtz operators. The first 
is simply as a discretization of ([!]), where the criterion is that the discrete solutions should 
approximate the true solutions well. Here we are particularly interested in the use of coarse 
meshes, say downto five or six points per wavelength, which are for example applied in 
exploration geophysics da eh ESI . The second application is internally in multigrid based 
solvers. In a multigrid method, the original mesh is coarsened by a factor two one or 
more times. On each of the new meshes a discretization of the operator is required. In 
this application the main criterion for a good discretization is that the multigrid method 
converges rapidly. The results concerning the application in multigrid methods are also 
of interest for recently developed two-grid or multigrid methods with inexact coarse level 
inverses which are currently some of the fastest solvers in the literature. (The 

method of [25J will actually be tested here.) Below we will write sometimes the fine level 
mesh for the original, uncoarsened mesh. 
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The small phase slowness errors for IOFD suggest that accurate solutions are possible 
even when quite coarse meshes are used, say downto five or six points per wavelength. 
We will show that this is indeed the case using numerical examples with constant, and 
smoothly varying velocity models (recall that k(x) = with c the medium velocity). 

We then consider the application of the IOFD discretization as coarse level discretization 
in multigrid based solvers. We will show that in this case IOFD can be used with very 
coarse meshes with downto three points per wavelength. With such meshes, solutions are 
generally not accurate enough for direct use, but the approximate solutions can still be 
used fruitfully in a multigrid method, where they are refined and iteratively improved. 
This is established using a set of two-dimensional examples, in which a two-grid method 
with IOFD at both levels converges rapidly (see also the results discussed in the next 
paragraph). As explained in [26], for the good convergence it is necessary to have very small 
phase slowness errors at these very coarse meshes. The IOFD method (in two and three 
dimensions) and the QS-FEM method (in two dimensions) are the only discretizations 
that have this property to our knowledge, and appear to be uniquely suitable for this 
application. 

In 3-D, the fact that a coarser mesh is used does not necessarily imply lower simulation 
cost. That depends also on the behavior of the solver. To investigate this aspect we 
present tests with a recently developed solver described in [25] . The solver uses a two-grid 
method with an inexact coarse level inverse, given by a double sweep domain decomposition 
preconditioner. As described in the previous paragraph, IOFD will also be used as coarse 
level discretization. Using the SEG-EAGE Salt Model with up to 10 8 degrees of freedom as 
example, we find that for downto six points per wavelength the cost per degree of freedom 
changes little when the frequency is increased. Computation time compare favorably to 
some of the results in the literature. 

The outline of this work is as follows. In section [2] the theory for finite difference 
discretizations of the Helmholtz equation with constant k is developed. The symbols and 
phase slownesses are defined and the discrete Green’s function is studied. In section [3] we 
consider the case of variable k and describe ray theory for discrete Helmholtz equations. In 
section [4] we compute the phase slowness errors of various existing schemes, as a reference 
for the new method. In section [5] we introduce our new interpolated optimized finite 
difference method. Section [6] contains some numerical simulations illustrating the accuracy 
of the solutions when using the IOFD discretization. Section [7] discusses the use of IOFD 
in multigrid based solvers. Finally, section [8] contains a brief discussion of some further 
aspects. 


2. Theory of discrete Helmholtz equations with constant k 
In this section we study finite difference discretizations of Helmholtz equation 
(2) Hu = /, H = -A - k 2 

in case k is constant. We will assume the grid is given by (h'L) d . In this and the next 
section it is convenient to write cq/3,... for multi-indices associated with grid points, such 
that with a = (oq,..., o^) is associated the grid point ha. A difference operator will 
be viewed as an operator on functions of x G ( hZ ) d . In this and the next section the 
dimension d can be any positive integer. 

For constant £;, a finite difference discretization of the Helmholtz operator H is a trans¬ 
lation invariant difference operator with coefficients depending on the grid spacing h and 
on k. By dimensional analysis we may assume that the matrix elements p a ^ of such a 
difference operator P are defined in terms of a finite set of functions fry by 

Pa,(3 = a—p{hk\ 


(3) 
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where / 7 is only nonzero for 7 in some finite set Sten(P). 

We will first consider the action of such an operator in the Fourier domain and define 
the associated symbol and phase slownesses. We next define a “dimensionally reduced” 
symbol. Then we consider the discrete Green’s function, i.e. solutions to the equation 

(4) Pu = 5, 

where the S function on (hZ) d is defined by 8{ha) = h~ d 8 ai: o ... 5 ad ^. We obtain the 
general solution of this equation in the Fourier domain and the asymptotics in the spatial 
domain, and determine the same information for the unique outgoing solutions. 

In the last part of this section we consider a modification of @ where first a function 
v is determined that satisfies 

(5) Pv = QS 
and then u is set equal to 

( 6 ) u = Qv. 

In this case we assume Q and Q are difference operators of order zero. Based on translation 
invariance and dimensional reduction, we assume that their matrix elements q a ^ and q a ^ 
are given by 

(7) q a , f 3 = g a -p(hk), q a ,p = g a -p{hk), 

where the g 7 , g 7 are smooth functions that are only nonzero for 7 in finite sets Sten(Q), 
Sten(Q). A solution u to such a system will be called a modified Green’s function. Our 
discretization of the Helmholtz equation will be a system of the form © and ([ 6 ]), where 
S is replaced by the right hand side /. 

2 . 1 . Symbol and phase slownesses. To define the symbol, we first define the forward 
and inverse Fourier transforms of a function u(x), x E ( hZ ) d . They are given by 

( 8 ) Fu{0 = h d Y, u(x)e~^ x 

X£(hZ) d 

(9) F~ l U(x) = ( 2 vr)“ d [ U{£)e* x d£, 

J [—7r/h,7r/h] d 

where the domain of Pu is [—7i/h,7i/h\ d . For constant k the finite difference operator P 
acts like a multiplication in the Fourier domain 

( 10 ) = 

where the function P(£), called the symbol , is given by 

( 11 ) p(o = h~ 2 YM hk ) eihl< - 

7 

This is similar to the continuous case, where the Helmholtz operator H acts by multipli¬ 
cation with H(£) = — k 2 in the Fourier domain. 

The Helmholtz equation has propagating plane wave solutions. These are functions 
u = e lx '^ that satisfy the homogeneous Helmholtz equation 

(12) He* x = 0. 

They are exactly the plane waves for which £ is in the zeroset Zh of the symbol of 
i7(£) = £ 2 — k 2 (This is of course the set of vectors of length k for H as defined, but the 
concept applies more generally.) If P is a translation invariant discretization of —A — k 2 
on M, d we can similarly look for vectors £ such that 

(13) Pe ix< = 0. 

These are the vectors in the zero set Zp of P(£). 
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If P is a discretization of the Helmholtz operator H then typically the set Zp is close 
to, but not identical to Zp. In other words, there are small differences in the wave vectors 
of the propagating waves, Zp ^ Zp. If Zp and Zp can be parameterized by angle, i.e. 

(14) Zp = {gp(e)6\eeS d - 1 }, 

and similar for Zp (for H(£) = £ 2 — k 2 this is of course the case), we define the relative 
wave number error as a function of 9 E S d_1 by 


(15) 




gp( 0 ) 

9h{Q) 


Closely related quantities are the phase slowness and phase velocity errors. If k = ^y, 
£ G Zp, there are associated phase slowness s p h and phase velocity vectors p p h given by 


(16) 


'Sph & 


Pph = 


IICII 2 


see e.g. [§]. The quantity S p h defined in (15) may hence also be called the relative phase 
slowness error, or simply phase slowness error. 

The actual phase error between a numerical and an exact solution is given by (see also 
subsection 2.3) 

(17) phase error = 2i— 

A 


where L is the distance between source and observation point, A is the wavelength and 5 p h 
is the phase slowness error associated with the particular angle of propagation. Because 
it is proportional to L/X the phase error easily may become dominant if is not careful in 
the choice of discretization in the high-frequency regime. 


2.2. Dimensional reduction and Helmholtz-like symbols. In case of coefficients, 
the symbol for arbitrary h can be expressed in terms of that for h— 1 

(18) P(0 = ^Pi(ht,hk) 

where 

(19) Pi(e,o = E e ^' 7 L( 0 - 

7 

By dimensional reduction the symbol Q(£), Q(£) can be written as 

(20) Q(0 = Qi(h£]hk(x)). 

where 

(21) Qi(70 = E e ^(0- 

7 

and similar for Q. 

To obtain the results below, we assume that the symbol P is Helmholtz-like as defined 
in the following 


Definition 1 . A symbol P(£) is said to be Helmholtz-like if the zero set Zp can be param¬ 
eterized as in (14), Zp is contained in ] — 7r/h,Tr/h[ d , P( 0) is negative, dP/d£ 7 ^ 0 at all 
points in Zp, and the map 


( 22 ) 


N :Z P ^ 5 d_1 : f ^ 


dp/om 

\\dp/dm\ 


that maps a point in Zp to the unit normal to the surface is a diffeomorphism. 


It follows that P is Helmholtz-like if P\ is Helmholtz like. 
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2.3. The discrete Green’s function. A Green’s function u(x) for the discrete equation 
will be defined as a solution of 0 - If k is constant the equivalent equation for the Fourier 
transform £/(£) reads 

(23) pm(o = i- 

We will first describe the general solution to this equation in the Fourier domain. Then 
we will consider the asymptotics in the spatial domain. Using the results obtained, we 
can then derive a unique outgoing Green’s function in the Fourier domain, and state its 
asymptotics. 


Due to the zeros of P, problem (23) has non-unique, distributional solutions. To explain 


their nature, we recall the closely related one dimensional problem to determine all / such 
that 

(24) xf( x) = 1. 

We also write this as M x f — 1, where M x is the multiplication operator by the function 


x. The solutions to (24) are the distributions of the form m 

(25) f{x) = p. v. - + b8, 

X 

where b is a free constant. Here the distribution p. v. ^ is defined by 

(26) 1 - — l 


(p.v. -,</>) = lim / 

x e ^° Jr\ 


dx. 


x 


In other words S is in the kernel of M x , while p. v. - is a particular solution to (24). 


In case of (23) we similarly have a nonzero kernel of Mp, with elements BSz P , where 


Sz P denotes the singular function of Zp, which is the distribution given by 


(27) 


Sz P {4>)= / 4>(x)dS(x), 


' Zp 


and B is any distribution on Zp. For functions / on W d with zero set Z such that 
V/(y) / 0 for all y £ Z. the principal value u = p. v. can be defined as follows. Let 

Z e — {x | d(x, Z) < e}, then 

4>{y) 


(28) 

We obtain 


(p.v. 


, 0 ) = lim / 

/ (y) e_s '° Jr 


R d \z e f{y) 


dy. 


Proposition 1. The solutions to fSS | ) are given by 

(29) U(0=P-v.jB + BS Zp 
where B can be any distribution on Zp. 

The freedom in the choice of B is related to the fact that in the spatial domain one can 
add any linear combination of plane waves e lx '^ with £ G Zp and still have a solution. 

Let u(x) be the inverse Fourier transform of a solution £/(£) for some smooth B 

(30) u(x) = (2ir)~ d [ p. v. ~^~e ix < d£ + (27r)~ d [ Be ix <. 

J [—7T /h,ir/h\ d -*Vs) J Zp 

We will study the asymptotic behavior of this integral for large ||x|| using the method 
of stationary phase. For p G we define a certain curvature-like quantity K(p) as 
follows. After rotating the coordinates, we may assume that dP/d^(p) is parallel to the 
d- th coordinate axis and that Zp is locally a graph 

(31) U = g(€i, ■ ■ ■ ,£d-i)- 
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By the assumptions g has a nondegenerate local maximum at (pi,... ,Pd-i). Let — Ai,..., —Xd-i 
denote the eigenvalues of the second derivative matrix d £ (pi,... ,pd_i). We de¬ 
fine 


(32) K(p) = \ 1 \ 2 ...\ d - 1 . 

For d — 3 this is the Gaussian curvature of the surface. In the following proposition TV -1 
denotes the inverse function of the map N defined in (22). The result and its proof have 
some similarity with results of Lighthill m • 


Proposition 2. Let u be the inverse Fourier transform of a distribution U as in 
such that B is a C°° function on Zp. Let = £ + (x) = N~ 1 (x/\\x\\) and £_ = £_(x) = 
N~ 1 (—x/\\x\\). The function u satisfies 


d+1 (d—l)7ri ^ ^ 

u(x) = 2 e—^\\x\\- — K^+y 




7T l 


Vi|5P/ae(f+)|| 


+ B{U) )e ix<+ 


(33) 


d+1 (d-l)ni 


+ (27r) 2 e 4 ||x 




Til 


\\dP/dm-)\\ 


+ £(£-) K" 


Jx-£- 


+ 0(\\x\\- 1 / 2 - d / 2 ), ||x||->oo. 


Proof. We start with the first integral in (30). For x G {hT) d the domain is really a torus 
and the integrand is C°° as a function on the torus. It is convenient to replace the integral 
on the torus by an integral over a bounded subset of M. d . Let t/r be a smooth, positive 
function supported in [—7r//i —p,7r//i + p], that is one on ] — Tv/h + g^n/h — g[ and satisfies 
YliZ-oo Vh( x + 27 rl/h) = 1 for x G M, and let 


(34) 

Then we can write 


(35) (2t r)~ d [ p. v. ^ ^)~ d 

J [—7v/h,7r/h\ d P\S/ 


ip( x ) = ■ ■■iPi( x d) 


l 


x G 


/ 

Jr* 


ip p. v. 


Ppe r(0 




di 6 


for x G ( hZ) d , where P per is the periodic extension of P, and this formula may also be 
considered for x G M. d . We assume g is sufficiently small such that Zp is supported in 
] - 7 r/h + t/, 7 v/h - g[ d . 

We will write x — tv, v G S d ~ l and consider the limit r -A oo. We assume coordinates 
are rotated such that v — (0 ,..., 0,1), using the same notation for the new coordinates as 
used so far for the old coordinates. 

The integral on the right hand side of (35) will be written as a sum of integrals over 
subsets using a partition of unity. For some smooth cutoff function y, denote 


(36) I x = ( 2 t r) d J x(£)U£) P- v - elTU 

We may assume there are four different types of x 

(i) x = x+ is one on a neighborhood of 

(ii) x = X- is one on a neighborhood of 

(ill) on supp x n Z P we can write Z P as a graph £ k = gfa,£k-i, €k+i, €d) 

(iv) supp x fl Zp = 0 

We consider these four cases in the limit r ^ oo using the method of stationary phase 
m- In case (iv) the integral I x = 0 (t n ) for any N by the lemma of non-stationary 
phase and we don’t need to consider this case further. In case (iii) we can write 


(37) 


X^p. v 


1 

' iWO 


= e (0 p- 


i 

6 : — s(£ i, • • • >£fc-i)£fc+ij • • • ,€d) 
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for some smooth function C(£) and perform the integral over This yields a smooth 
function of (£ 1 ,..., £&_i, £&+!, €d)- By the lemma of non-stationary phase it follows that 
again I x = 0(r~ N ). 

In case (i) we can write Zp locally as a graph ^ i,..., £d_i). For brevity denote 
— (£ 1 ,..., £d-1 ) We observe that we can write 

(38) xV’P-v. -p = ho(0 + hi(g) p. v. — - MUK 1 - ~ g(O)) ^ ^ • 

where are smooth, compactly supported functions and ^2 = 1 around 0 and 

/c(£+) = ||ap/ag(C+)|| • Then f° r the t erm the lemma of non-stationary phase can be 

invoked. Hence this term is 0{r ~ N ) for any N. For the third term, the same result can 
be obtained using integration by parts. For the second part we recall the standard Fourier 
transform ^p. v. = = —i7rsgn(^), it follows that 

(39) =e ia y l -sgn(y). 

rj — a Z 

As a consequence, we obtain 

(40) I x+ = (27r) - ^ -1 R J MOe iTfl( ° < + 0(t~ n ) 

any TV. We can now apply the stationary phase lemma. The function g has its maximum 
at and can be expanded as, possibly after a further rotation of coordinates 


d -1 


(41) 


g(£i, ■ ■ -:U- 1) = v ■ £+ ~ T VU _ (Oi) 2 + o(ll£' - A 

i=i 


see the discussion preceding (32). This yields 
(42) I = j L(2 7 r)-( d - 1 )/ 2 _ - _ 

( j x+ 2 ( } ||ap/ae(c+)ll 

The contribution I x _ in case (ii) can be computed similarly, resulting in 


e -(d-l)^ e iTV-i +T -{d-l)/2 K ^ + yl/2_ 


(43) 


4- = -o( 2 ^: 


r(d-1)/2. 


1 


\\dP/d^-)\\ 


e (d-l)^ e iTt)-{- r -(d-l)/2_ K -^_^-l/2_ 


For the surface integral 
(44) 


(2ir)~ d [ B(d)e ix< dS(ti), 

J Zp 


we again assume x — tv and consider the limit r -A oo. A partition of unity is applied and, 
by the method of stationary phase, the only contributions that are not 0{r~ N ) for any 
N come from neighborhoods of £±. To determine the contribution from a neighborhood 
of £+, we assume that u is parallel to the d-th coordinate axis so that Zp is locally given 
by a graph ^ = #(£'), = (£ 1 ,..., £d-i)- The method of stationary phase can be applied 

directly. The only contributions come from neighborhoods of £±, and can be computed 
similarly as for the integral 


The contribution 4 + , 4- and the two contributions from the integral (H4| together 
give the result. □ 

It is straightforward to obtain the outgoing solutions to Q and (23). In (33) the term 
with phase factor e %TV '^~ must vanish, and we obtain the equation 


(45) 


m 


iri 

\\dp/dmr 


for f G Z P . 


We state this as a theorem and include the asymptotic expression for the solution in the 
result. 
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Theorem 3. The outgoing solution to (23) is given by 


(46) 


U(0 = P-v. 


7T iS Zp (0 

p(o ' \\dp/dm\' 


+ 


d-1 (d-l)iri d-1 iK(£-1 _) 


2 e 


x 


II^M(U)II 


e ix <+ + 0(\\x\\-P 2 - d / 2 ), 


Its inverse Fourier transform u(x) satisfies 

(47) u[x) = (27r) 
where K and are as in proposition [S| 

The above analysis can be repeated for continuous, Helmholtz like operators with the 
same result (see also m)- For the usual Helmholtz operator H in d — 3 dimensions we 
have that Z H given by ||£|| = k, \\dH/d£\\ \^ e z H = 2 k, K + = ^ and 

(48) „( I ) = _J_ e «ll*ll + 0 (||x||- 2 ), 

so that the highest order asymptotic expansion actually equals the well known outgoing 
Green’s function. 

2.4. The modified discrete Green’s function. Let 

(49) H 1 (t,k) = e~k 2 . 

It follows from theorem [3] that if Hi and Pi have the same zero sets, i.e. identical phase 
slownesses, then the solutions to Hu — S and Pu = S have asymptotically the same phase. 
The amplitudes however will differ by a factor evaluated at the zero set. In this 

subsection we consider therefore the solutions u to the equations ^ and ([ 6 ]), which, as 
we will see, obtain different amplitudes. 

The Fourier transformed solution £/(£) to ([BJ and © is given by the product of the 
solution U given in proposition [l] and a factor Q(£)Q(0- Using this, we can formulate a 
result similar to Theorem [3j In this case the adjective outgoing refers to the solution v of 
©• The result can be proven by similar arguments as used to prove proposition [ 2 ] and 
theorem 03 

Theorem 4. The Fourier transform of the outgoing solution to m and & is given by 


(50) 


a(0 = Q(0Q(0 p-v. 


+ 


niSzp(£) 


V 'go nwyaecoir 


Its inverse Fourier transform u(x) satisfies 
(51) 

" " \\dP/d^ + )\\ 


e ix <+ + 0{\\x\\-P 2 - d / 2 ), 


u(x) — (27 t)~ 
where K and are as in proposition [S| 

Summarizing our findings so far, the discrete solutions u to © and © are asymptoti¬ 
cally equal to the solutions of the continuous Helmholtz equation Hu = S if the following 
two conditions are satisfied 

(i) Pi(£, k) and FZi(£, k) have the same zero sets 

(ii) Q i and Qi satisfy 

(52) Q,(( k)OM k)~ 

(52) || OTl/aeKjS .)|| 

for all (£, k) such that Pi(£, k) — F/i(£, k) — 0 
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3. Theory of discrete Helmholtz equations with variable k 


In this section we define a class of discrete approximations to the Helmholtz operator 


with variable fc, together with the associated symbols. This is the topic of subsection |3.1 
We then study ray-theoretic solutions to the equation ® and to the set of equations 
0 > where P and Q are now variable coefficient operators. 


We assume that k = -tA, where c is smooth and we consider the limit u 

c{x) ’ 


discrete case we assume that huo — constant, 
the ansatz 


oo. In the 

Ray-theoretic solutions are then based on 


(53) 


u(x,u) = A(x, 


for some smoothly varying A and 4>. For the continuous Helmholtz equation, such solutions 
are well known and are constructed in two steps. First the ansatz (53) is inserted in 


the PDE, and an expansion in cj is performed. Requiring that the highest order terms 
vanish leads to the eikonal equation for and the transport equation for A. Secondly, 
initial/boundary conditions for these equations are obtained from the asymptotic behavior 
of the constant coefficient solutions. In this way, the solution modulo an error of lower 
order in uo is obtained. 

For our class of difference equations we follow the same program. The constant co¬ 
efficient solutions were already analyzed in subsection |2.3[ In subsection |3.2| we find a 


nonlinear first order PDE for $ and a transport equation for A. Remarkably, we obtain 
the same equations in the continuous and discrete case when formulated in terms of the 
symbols (which are defined for both continuous and discrete problems). See [19] and m 
for the continuous case and methods used in that case as well as here. 

In the last part of this section we consider the ray-theoretic solutions to 0 and 0 - 
The conditions ^ and ([h]) from subsection 2.4 for P and Q to obtain accurate solutions, 
need to be modified and extended to have the same ray-theoretic phase and amplitude in 
the continuous and discrete case. The operator P should be discretized using a symmetric 
discretization (with = 1 / 2 , see below) and we should have Q = Q. This is the topic of 
subsection 13.31 


3.1. Symbols and operators for variable k. In case k depends on x, finite difference 
discretizations of the Helmholtz operator may depend in different ways on the function 
k. For example, the coefficients p a ^ may depend on k and its derivatives at x = ha , but 
they may also depend on k at different points, for example on k{ha) and k(h/3). We will 
consider a class of difference operators P, where the matrix elements p a ^ depend only on 
the value of k at (1 — t)ha + th(3, where t E {0,1/2,1} is a fixed constant. In other words 
we consider operators P with matrix elements of the form 

(54) p a ,p = jpf a _p(hk(( 1 - t)ah + tj3h)). 

Note that the operator is symmetric if t = 1/2 and / 7 = /_ 7 . This will turn out to be an 
appropriate choice for a discrete Helmholtz operator. We will assume that k(x) is defined 
for all x, not only those in the grid. Similar we assume that for Q we have 

(55) = g a -p(hk(( 1 - t)ah + t/3h )). 
and similar for Q. 

For such operators it is not obvious how to define the symbol. To find an appropriate 
definition, we first consider how to define an operator from a symbol H(x,£) in the con¬ 
tinuous case. This is the subject of pseudodifferential operator theory, and can be done 
with the formula cun 

(56) Op t (H(x,£))u = (2n)~ d JJ H(x + t(y - x), £)e* (a;-J/) 'U(y) dy. 
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A map from a function H(x,£) to an operator such as Op t {H(x,£)) is called a quanti¬ 
zation. For t — 0, the previous formula is the standard left-quantization, t — 1 is the 
right quantization and t = 1/2 is the Weyl quantization. If = £ 2 — fc(x) 2 , then 

Op t (H(x,£)) = —A — fc(x) 2 , independently of which of these quantizations is used. 

To obtain a symbol associated with the operator P defined in ( |54| ) we rewrite the 
expression for Pu(x) as follows 


(57) 


Pu(x) = h 2 f^(hk(x + t r yh))u{x + hj) 

7 

= h~ 2+d Y^ fj(hk(x + t(y - x)))S(x + hj -y)u(y). 

7 ye(hZ) d 


Using the Fourier domain representation 5(x) — (27r) d S[-n/h n/h\ d e%X ^ this can rewr ^ _ 

ten as 


(58) Pu(x) = h~ 2+d { 27r)~ d V [ y2Mhk(x + t(y-x)))e^ x+h ^<d^. 

ye(hZ)d d[—Tr/h,Tr/h] d 7 

This can be written in similar form as ( |56| ), namely as 

(59) Op t (P(x, £))u(x) d = (2vr )~ d V [ P{x + t(y - x),^)e^ x ~ y)< u(y) d£, 

ye(hZ)d J[-n/h,n/h\* 


where 


(60) 


P(x,0 = h 2 '£2f 1 (hk(x))e lhl< . 


Thus, associated with P defined in (54) is associated the symbol P(x,$) given in (60). 
The parameter t corresponds to the type of quantization, left, right or Weyl quantization. 

With these definitions, the symbol (60) for variable coefficients may also be expressed 
entirely in terms of P\ 


(61) 


P(x,Z) = ±P 1 {ht;hk(x)). 


A symbol P(x,£) is called Helmholtz like if it satisfies the definition for each fixed x. 


3.2. Ray-theoretic equations for amplitude and phase. In this section we consider 
the high-frequency limit uo -A oo. We assume that coh = constant, and recall that k(x) = 
where c(x) is C°°. The operator P and the symbol P(x,£) become cj-dependent. By 

P(x,£) we denote the symbol for uo — 1. 

(62) 

For other values of uo we find that 

(63) P(x, a;) = uj 2 P(x , —) 

(jj 

We consider the action of P on functions of the form 

(64) u{x) = e iu ^A(x) 

where and A are C°° functions. From the symbol P(rr, £;cj) and the phase function 
one can derive naturally a vector field, which we call (cf. [TO, section 4.3]) 

dP 

(■ L P ,$,«;). = —(a:,wV$;w). 

J °Zj 


( 65 ) 
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This vector field is determined by Lp ^ = V<h) as follows 

( 66 ) L P ^ tU) (x) = uLp^{x). 

Proposition 5. We have 

(67) 

e -i U 9(x)p^(x) A ^ v$(x))A(x) 


1 f \\. . dA 1 . . . . . ^^ d 2 P 

+ oT + ~( dlvL P,$) A +(*- !/ 2 ) 52 


( 68 ) 


$7 + y) = $(z) + V$(z) • y + i dx^x k ViVk + 


i \ ^ ^ J dxj 2 

J 3 

T 0(1), —)• oo. 

Proof. The proof uses a Taylor expansion of the phase function to second order 

i^2 92$ 

j,k 

a Taylor expansion of the amplitude to first order 

(69) A(x + y) = A{x) + Vi(s) • y + 0(||y|| 2 )) 
and a Taylor expansion of the matrix coefficients to first order 

(70) p a ,p = ^ fa-p(hk(ha )) + ^ f' a _p(hk(ha))hVk(ha ) • h(/3 - a). 

The exponent e ia,$ ( x +s/) ig then written as a product of three factors 

(71) = e iw#(*) e iwV*(*)-y j 1 + f^rViVk + OMMI 3 ) 

\ 3 

These expansions are inserted in the sum 

(72) (Pu)(x) = Jf 2 f'yihkfa + th^))u{x + / 17 ). 

7 

The factor e za;< ^) can be put in front of the expression outside the summation 


dx^ 


(. Pu)(x ) =e^W — 


2 f^hk(x))A{x) 

7 

+ (£ e iwv*(*)-fc7 p{hk{x))yA{x) ■ (hj) 


(73) 


+ liwA 


<9 2 $ 


—' dxjdxb 

j,k 


(hpj{h^)kfj{hk(x)) 


+ A te iujV ^ x) ^f^hVk ■ (hj) 


+ 0(h 2 ) 
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We next use the expression for the symbol P(x,£) = ^4 X 7 f-y(hk(x)), the following 
expressions for the derivatives of P(x,£) 


(74) 

(75) 

(76) 


dP _ i 
d 2 P 


f y (hk(x)) 


dxjd^k h 2 


= P k ( x )) h 


dk 

8x; 


d 2 P 

i-KM k 


~ 5^(^7)i(^7)fce* fcr{ / 7 (/ifc(x)). 


and an expression for the derivatives of 
d(Lp^)j _ d 2 P 


(77) 


This yields 


dx k 


dx k dP 


(x, wV<h) + u 


d 2 P d 2 § 
d£jd£i dx t dx k ' 


(Pu)(x) =e iuj *(*> 


(78) 


P(x, wV$(s);w)4(i 

i / n r)A 1 1 r) 2 P 

+ 7 ( + 2 div L Pp.w A + (* - h) 

\ j —1 3 

+ 0 ( 1 ) 


2 , d x j€j 


UJ —)► OO. 


Using equations (63) and (66) the result follows. 


□ 


The result is similar to the result in Proposition 4.3.2 of m- To find A and 4> such 
that 

(79) P(e iu ^+A(x)) « 0 

the phase function 4> must satisfy the equation 


(80) 


P(x, V4>) = 0, 


which is a nonlinear first order equation like an eikonal equation, and the amplitude must 
satisfy a transport equation 


(81) 


O A 1 f) z P 

: + x(div Lp^)A + (t- 1/2) y] n at = 0 


a 2 p 




r' fajdtj 


For t = 1/2, this equation conserves |yl| 2 . 

Ray theoretic solutions to equation 0 can now be constructed just as in the continuous 
case. By a rescaling, theorem [3] can be used to obtain the asymptotics of a solution for 
x ^ 0 and uo oo. The amplitude and phase from formula (47) can hence be used as 


initial/boundary values for the eikonal equation for 4> and the transport equation for A, 
and these 4> and A can be determined from these equations, where we note that the eikonal 
equation may not have globally defined solutions, just as in the continuous case. 

We briefly recall the continuous equivalent of Proposition [5j The following is basically a 
reformulation of proposition 4.3.2 of m and can be proven using the method of stationary 
phase found in the same text. 
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Proposition 6 . Let H be a continuous Helmholtz like symbol. For the action of Op (H) 
on e luJ ®A(x) we have the asymptotic development 
(82) 

e -iu<s>{x) Op t (H)(e l ^ x) A(x)) =lu 2 H(x,V$(x))A(x) 

1 /v^ /r . dA 1 d 2 H 

+ ^ + 2 dlv L H,*) A +{t ~ V 2 ) 5^ d^ A 


+ 0 ( 1 ) 


UO 


OO. 


3.3. Amplitude correction. In this section we consider ray-theoretic approximations 

v = e luJ ^ B and u = e ZUJ ®A for the solutions v and u to ^ and @> . Assume we have 

reference ray-theoretic solutions u Te f = e tu} ^ ref ^ A re f(x) associated with Helmholtz equa- 

2 

tion Hu re f = 5 where H(u) = —A- rx 2 -> i- e - ^ref satisfies the eikonal equation and 

C\X) 

A Te f the transport equation with appropriate initial conditions. In the following we let 
Hi& k)=e~ k 2 , H(x, 0 = ( 2 - 
Assume that 

(i) Pi(£, k) and i2i(£, k ) have the same zero sets 

(ii) Q i and Q\ are identical and Qi(£, fc) = Qi(£,k) = Qi(£, fc) satisfies 

Wdp^oaem 


(83) 


Qi(e,fer = 


woih/d^m 

for all (£, fc) such that Pi(£, k ) = fc) = 0; 

(iii) P and H are derived from their respective symbols using t — 1/2 quantization . 

We argue that in this case, to highest order u has the same ray-theoretic approximation 
as u re f. We omit a formal proof, because the arguments are similar as those used above. 

The construction of the phase and amplitude functions T and B proceeds almost in 
the same way as for solutions to ©• Eikonal and transport equations are as follows from 
Proposition |5j The constant coefficient solutions differ by a factor Q(£+) and have the 
same phase, resulting in different initial/boundary conditions, such that on a small sphere 
T around 0, where we impose the initial/boundary conditions for the eikonal and transport 
equations, we have 

. . ^(*) = $ref(*) 

B(x) = Q u= i(x, V$ re f) Tef(z) 

As a result, 4/(x) = 4> re f(x) everywhere. While we have different transport equations the 
operators Lp ^ f and Lp^ f are scaled versions of each other 

(85) Lp ^ ref = Qw=i(x, V<h re f) P#- 5 $ ref5 i 

It follows from this fact and the transport equation for A re f, that 

9QZh A rct 


( 86 ) 


Eu 


+ (divL^JQj^ = 


^ref 


dxj 


= 0 . 


This and (84) shows that 

(87) B(x) = Q~ l (x , V$ re f) A lef (x) 

everywhere. 

The function u is given by applying Q to v. The action of Q on e luJ ^B is to highest 
order equal to a multiplication by Q(x, cj), so that 

= 4/ = 4> re f and 
A{x) = Q(x, V$ ie f)B(x) = A re{ (x), 


( 88 ) 
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concluding the argument. 

4. Phase slowness errors for existing discretizations 

In this section we will describe three types of discretizations of the Helmholtz equa¬ 
tion 0. namely standard finite differences, compact finite differences and Lagrange finite 
elements on regular meshes. We then compute phase slowness errors to compare the per¬ 
formance of the different methods in this respect, and to obtain reference values for our 
new method constructed below. 

We modify the notation compared to the previous two section. In this section the degrees 
of freedom for all three types of methods are denoted by (in three dimensions) and 
associated with a regular mesh with grid spacing h. For finite element methods of order N 
the cells are of size Nh, and cell boundaries are located at j, fc, l = 0 mod N. Occasionally 
we will use d — 2 or 3 to denote the dimension of space. 

4.1. Standard finite differences. In a standard finite difference discretization of the 
operator —A — k 2 each of the one-dimensional second derivatives in the Laplacian A = 
da? + fa? + fa? * s a PP rox i ma t e d by a central difference approximation of the given order. 
These are given by 

N/2 

(89) D { 2 N) ui = h~ 2 ^ cWui +m 

m=—N/2 

where the c are as in the following table for N = 2, 4, 6 , 8 



m — — 4 

-3 

-2 

-i 

0 

i 

2 

3 

4 

05 

II 




1 

-2 

1 




4 



1 

4 

5 

4 

1 





12 

3 

2 

3 

12 



6 


1 

3 

3 

49 

3 

3 

1 



90 

20 

2 

18 

2 

20 

90 


8 

l 

8 

1 

8 

205 

8 

1 

8 

1 

560 

315 

5 

5 

72 

5 

5 

315 

560 


The discrete approximation to the term —k(x) 2 u in (JTJ) is simply given by —kf rnn u^ mn . 
The two-dimensional case can be done similarly. 

4.2. Compact finite difference discretizations. For constant P, compact finite differ¬ 
ence discretizations take the form 

(90) (Au)i ,m,n — ^ ^ ^p,q,r u l-\-p,m-\-q,n-\-r • 

(p,q,r)e{- i,o,i } 3 

Because of symmetry, there are four different coefficients Aj, j — 0,1, 2, 3 and 

(91) a p,qx ^IpI+M+M? (p* 0 ^ { — 0,1} 

In 2-D we have 

(92) {Au) t ,m ^ ^ Q'p,q'U'l-\-p,m-\-q m 

(P, 9 )€{- 1 , 0 ,l } 2 

and there are three different coefficients A/, j = 0,1, 2 and 

(93) a Ptq = A \ p | + |,|, (p,q)e{- 1,0, l} 2 . 

The choice of coefficients is done in different ways in [3} [13i QJU 21, 27, [6, 29, [26]. The 
QS-FEM method [3] is a two-dimensional method, for which the coefficients are given 
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modulo an overall normalization by 
^QS-FEM = 4 

4 QS-fem _ 0 ci(a)si(a) - 02 ( 01 ) 52 ( 01 ) 

(94) 

C2(a)s2(a)(ci(a) + si(a)) - ci(a) 5 i(a)(c 2 (a) + 52(a)) 

^qs-fem __ c 2 (a) + 52 (a) — ci (a) — si (a) _ 

2 C2(a)s 2 (a)(ci(a) + si(a)) - ci(a)si(a)(c2(a) + 52(a)) 

where a = kh and the auxiliary functions ci, si, C 2 , S 2 are dehned by 


(95) 


ci (a) 
c 2 (a) 




cos 



We will not discuss the fourth order method of H3J because it contains still a free parameter 
and one of the authors has later published a sixth order method in m- In the latter 
method variations of k are taken into account. In case of constant £;, the coefficients are 
given in three dimensions by 


(96) 


^CHOe = + M _ 
0 15 

/1CH06 _ 7 

Al - ~ 15 + 


14 k 2 h 2 k 4 h 4 

+ 


and in two dimensions by 


(97) 


4 CH 06 

^0 


/1CHO6 _ 

" “ 6 ' 


15 

k 2 h 2 

TcT 


10 

y 


20 

a cho 6 


1 

10 


2 r.2 


k 2 h 

”90” 


41 k 2 h 2 k 4 h 4 

+ 


45 

k 2 h 2 

~90~ 


20 

^CKoe 


/1CHO6 


k 2 h 2 
~90" 


1 

30 


again modulo an overall constant. For the method of Sutmann m we have (this method 
is only for 3-D) 


4 UT = (1 - \ kV + m kihi ~ tA.^ 6 ) 

(98) 15 

4 UT = - 7 (1 - iiV) 4 UT = (1 + ±kV) 


4 sut 

A 3 


1 

30 


In [16, | 2 Tj 6 . 26] the contributions to —A — k 2 are split in a contribution from —A and 
a contribution from — k 2 . We define a three dimensional, symmetric discretization of the 
identity M depending on three parameters ay, oq, a ^3 by 


(Mu)i^ n / 'HIp,g,r^'/+p,m+g,n+r 

(99) (p,^,r)e{-i,o,i} 3 

m p , q , r = M lp]+lq]+lr] , (p, q, r) G {-1,0, l } 3 

where now Mq = ay, M\ — M 2 — f§, M 3 i-exi-a 2 -a 3 . p or cons tant k the discrete 
form of the term —k 2 u is given by —k 2 {Mu)i^ m ^ n . Before defining the negative Laplacian 
we define a two-dimensional weighting operator, discretizing the identity, depending on 
two additional parameters 014,015 


(7V^ ’ ^)/ ?m?n ^ ^ 

p,q£{- 1 , 0 ,l} 2 

™p,q N\p\+\q\> (.P? (?) 9 , 1 } 


( 100 ) 
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with TVq = 04 , Ni = N 2 = 1 ~ a ^~ as . Each of the second deratives in the Helmholtz 
operator will be discretized using the tensor product of a two-dimensional weighting oper¬ 
ator that discretizes the identity and the standard second order discrete second derivative 
D 2, where j — 1,2 or 3 indicates along which axis the second derivative operators. The 
resulting matrix is 

(101) - d[ 1] 8 > at! 2 - 3 ] - D® ® N l 1 ’ 3 ] - £>Jj 3] 8 iVl 1 ’ 2 ! - k 2 M 

In the two-dimensional case there are in total three parameters or, D 2 , <^ 3 , with Mq = aq, 
Mi = M 2 = 1 ~ a i~ a2 ? and TVo = « 3 , Aff = The coefficients Aj for the 3-D case 

are given in terms of the aj (modulo an overall constant h~ 2 ) by 

( 102 ) 

Aq = 6 D 4 — (kh) 2 a i A 2 = — \ct$ + ^(1 — oq — 0 : 5 ) — (fcft) 2 ^a 3 

A\— - aq + - ( kh) 2 \a2 A 3 = - |(1 — 0*4 — a 5 ) - {kh) 2 \{l - ol\ - ol 2 ~ a 3 ). 

For the 2-D case we have 
(103) 

Aq 4(^3 — (/c/l) 2 Di Aq = 1 — 2 d 3 — (fc/l) 2 |(T2 ^2 = — 1 + O3 — (fc/l) 2 |(l — Oq ~ OL2). 

An advantage of this formulation using tensor products is that in case of PML layers 
aligned with the coordinate axes, the second order operator can simply be replaced 
by its PML-modified versioiQ 

In this way we have derived a family of second order accurate discretizations. In m 
and | 21 | the same family of discretizations is considered in two resp. three dimensions 
(but differently parameterized), and coefficients are chosen such that the maximum phase 
slowness error is minimized, where the maximum is taken over all angles and a range of 
kh corresponding to at least four points per wavelength. This leads to the choices 

a O pT = 0.4964958 o£ PT = 0.4510125 ap PT = 0.052487 

'''' 1 ' a2 PT = 0.648355362 a^ PT = 0.296692332 

for the method of m and 

(105) af ss = 0.6248 af s = 0.37524 af s = 0.77305 

for the method of m ■ 

In [ 6 ] and [26j it is observed that smaller phase slowness errors are obtained when the 
parameters aj are allowed to vary. In [ 6 ] a set of 7 parameters (in three dimensions) 
is chosen piecewise constant. We will not describe this method in detail but refer to the 
paper for resulting phase errors. In [26j the above described set of 5 parameters are chosen 
as piecewise linear functions. However, in this work, the aim is different, because the phase 
slowness differences with a fine scale operator are minimized, not with the exact operator, 
so that the values of the phase slowness errors cannot be compared. 

4.3. Lagrange finite elements on regular meshes. For the description of Lagrange 
finite elements on regular meshes, which will also be used in some of the numerical exam¬ 
ples, we start with the one-dimensional case. In this case the finite element cells are the 
intervals ((j — 1 )Nh,jNh), j = 1 , 2 ,..., each containing N — 1 interior points and two 
boundary points. The reference cell is (0,7V), and shape functions on this reference cell 

bn a PML layer, say a layer associated with x\ — constant, the derivative is replaced by 
a <apML,i(^i) afy where <apML,i(^i) — i+z C r 1 (a;i)/aI an( l function g\ indicates the local amount of 
damping H3 H 0 . We choose 04 quadratically increasing. The discrete second derivative in the 
first coordinate in this PML layer becomes h _2 apML,i(xz, m , n )(<apML,i(^+i/2,j,A ; )(^z+i,m,n - ui iT n,n) ~ 
apML,i(^-i/ 2 ,j,fc)(^,m,n — ui -By rescaling the equations with a factor c^ml x the symmetry of 
the system is restored. 
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are given by standard Lagrange polynomials, which we will denote by L^ N \x), and which 
are one at x — j and zero at x — 0,1,..., j — l,j + l,...N — l,N and defined to be 0 
outside [0, N}. Letting fcGZ and l E {0,1,..., N — 1}, the one-dimensional trial and test 
functions are given by 


(106) 


^kN+l( x ) 


4 Ar) (|-(fc-l)AT) + L^ ) (|-fciV) if 1 = 0 
l[ N \% — kNh ) otherwise 


The two and three-dimensional trial and testfunctions are given by tensor products of the 
i/jj N \ The finite element discretizion is of course derived from the weak form 


(107) 


$(r 


n d 

,») d ='/E 


dudv f 2 7 f r 1 

— —— —ax — / k uv ax — / jv ax. 
oxj oxj Jq Jq 


The elements of the matrix in the discretization are given by 


( 108 ) a Z,m,n;p,^,r — ^{^Pp,q,r 5 V^,m,n) 

The contribution from the term J2j=i dx can be called the stiffness matrix and 

the contribution from J ^ k 2 uv dx can be called the mass matrix. If k is constant (or cellwise 
constant) , the stiffness and mass matrices can be computed exactly. If k is variable, then 
only the stiffness matrix can be computed exactly, and for the mass matrix some sort of 
quadrature must be used. For constant k these computations are standard and easily done 
using a computer algebra system, and we will not write down the resulting coefficients. 

For constant fc, the finite difference methods are obviously translationally symmetric, 
i.e. if we denote by a^ m n . p q r the matrix elements we have 

(109) a t 

,m,n]p,q : r — ( dl+A,m-\-B,n+C\'p+A,q-\-B,r+C 

For the finite elements there is a symmetry under a subset of translations given by the 
A, B,C that are multiples of N. 


4.4. Phase slowness errors. For finite difference methods, finding the phase velocities 
or slownesses comes down to determining the zeros of the symbol P(£) associated with 
a difference operator P, see subsection |2.1| The symbol is not difficult to obtain, for 
example, for compact finite difference discretizations, the symbol is 
( 110 ) 

P(£) = h 2 [A 0 + 2^4i(cos(/j£i) + cos(/ i£ 2 ) + cos(fo£ 3 )) + 2^4 2 (cos(/j(£i + £ 2 )) + cos(/i(£i - £ 2 )) 
+ cos(/i(£i + £3)) + cos(/t(£i - £3)) + cos (h(& + &)) + cos (h(& ~ 6))) 

+ 2^ 3 (cos(/i(^i + & + 6)) + cos(/i(^i - 6 + 6)) + cos(/i(^i + £ 2 - 6)) + cos(/i(^i - & ~ 6)))] 

To compute the zeros numerically, the standard numerical solver fsolve from Matlab was 
used, as well as the more accurate version vpasolve. 

For finite element methods the elements of the kernel of the operator are no longer 
simple plane waves, but Bloch waves. In the appendix is described how we compute the 
phase slowness errors in this case. 

Phase slowness errors are directionally dependent, i.e. they depend on 9 G S d_1 as 
explained in section [2T| They also depend on kh or equivalently on the number of points 
per wavelength G = We have computed the maximum relative phase slowness errors 
over 6 G S d ~ l for a number of schemes as a function of 1 /G. In two and three dimensions 
these schemes are the finite element schemes of order 1, 2, 3,4, 6 and 8 , the standard finite 
difference discretizations of order 2,4, 6 and 8 and the sixth order compact method of [29] . 

In the graphs below these results will be indicated by the letter FE1, FE2 etc., FD2, FD4, 
etc. and CH06. In two dimensions we also included results for the QS-FEM method of 
0 and the method of Jo, Shin and Suh [Kj, denoted by JSS. In three dimensions we also 



A DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 


19 




-FE1 


-FE2 


- FE3 


FE4 


- FE6 


- FE8 

-B— 

-FD2 

-B— 

-FD4 

-B— 

-FD6 


FD8 


■ JSS 

_ 

■ CH06 


n.Q-PPM 



Figure 1. Phase slowness errors for some 2-D schemes as a function of 
the inverse number of points per wavelength 1/G. 
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Figure 2. Phase slowness errors for some 3-D schemes as a function of 
the inverse number of points per wavelength 1/G. 


have the method of Operto Et A1 m , indicated by OPT4 and the method of Sutmann 
mi indicated by SUT. We have not included results on the method of [ 6 ], these are given 
in Figure 2(c) in that work. The phase slowness errors as a function of 1/G are plotted in 
Figure |T] and [2j At the end of section [5] we will briefly discuss these results. 

5. A DISPERSION MINIMIZING SCHEME WITH AMPLITUDE CORRECTIONS 

In this section we will define our new discretization of the Helmholtz equation. In this 
scheme, the approximate solution to the Helmholtz equation Hu — f , H — —A — k(x ) 2 is 
found by solving a discrete system 

( 111 ) Pv = Qf 
and then setting 

( 112 ) u = Qp, 

where P and Q are compact finite difference operators defined momentarily. 
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In section [5] we studied ray-theoretic solutions to difference equations of the type (111) 
and ( 112 ) with / = 5, and we observed that the ray-theoretic solution to these equation 


would be identical to those of the Helmholtz equation 
(113) Hu = 6, H = -A - k(x) 2 


if the requirements (i) to (iii) of subsection 3.3 are satisfied. It follows from the derivations 
that if these properties are not satisfied exactly, but there are small differences between 
the zero set of P\ and that of H\ and between the values of Qi(£, k ) 2 and then 

there will be small errors in the phase and amplitude of the ray-theoretic solutions. The 
operators P and Q will be chosen such that these differences are minimal. We will first 
construct P in subsection |5.1| Then Q will be constructed in subsection|5.2| In subsection 


we will discuss the phase errors of the new method. 


5.1. IOFD discretization of the Helmholtz operator. In sections [2] and [3] a general 
form for P was given in terms of functions f 1 of fc/i, see Q, ( |54| ). For the 5 or 3 dimensional 
operator family of subsection |4.2| (in 3 and 2 dimensions respectively), these are given by 


(114) 


{ 604 — (kh) 2 a\ for 7 G { — 1, 0, l} 3 , |y| = 0 

-04 + as - ( kh) 2 \a 2 for 7 G {-1, 0, l} 3 , jyj = 1 

-|o 5 + 7^(1 - 04 - 05 ) - (kh) 2 j^a 3 for 7 G {-1, 0, l} 3 , jyj = 2 

-|(1 - 04 - < 25 ) - {kh) 2 \(l -a 1 -a 2 ~ 03) for 7 G {-1, 0, l} 3 , jyj = 3, 


in 3-D and by 


( 403 - (kh) 2 ai for 7 G {-1,0, l} 2 , I 7 I = 0 

(115) / 7 = < 1 - 2o 3 - {kh) 2 \a 2 for 7 G {-1,0, l} 2 , jyj = 1 

[ -1 + 03 - {kh) 2 \{l - 01 - o 2 ) for 7 G {-1,0, l} 2 , | 7 | = 2. 


in 2-D. Here I 7 I = 1 71 1 +... +1 7 ^| and we used equations ( 102| ) and ( |103 ). We let aj depend 


on = 1 /G, where G is the number of points per wavelength used in the discretization 


(116) o J =o J (l/G), l/G = —. 

Next we will choose a parameterization for these function and we will describe how, by 
minimizing the phase slowness errors in a least-squares sense, we obtain suitable choices 
of the functions o^, j = 1 ,..., 2d — 1 . 

In [26j the aj where chosen to depend piecewise linearly on 1/G. Here we let aj depend 
piecewise polynomially on 1/G, using Hermite interpolation. We will specify a number 
of control nodes, and at each node the value of aj and its first derivative d(i/ J G) are 
prescribed. We will assume that the coefficients ay vary slowly, so that we can indeed 
define the four coefficient of the stencil using five parameters depending on 1/G. If nc 
denotes the number of control nodes, in this way the functions aj are parameterized by 
2nc parameters. In total we have (4d — 2)nc parameters, collectively denoted by P. 

Next we specify the objective functional. The first contribution to the objective func¬ 
tional is the square integrated phase slowness error, integrated over angle and 1 /G. 
Because of the symmetries, the phase slowness error need not be integrated over all 
e € s d ~\ but can be integrated over a subset of the sphere. In 2 dimensions, the 
angle variable can be chosen in ©2 = [0, 7 t/ 4 ]. In 3 dimensions, using spherical coordinates 
0 — (0i,02) — ( polar angle, azimuthal angle), the domain is ©3 = [ 0 , 7 r/ 2 ] x [0, tt/ 4]. The 
second contribution to the objective functional is a regularization term involving ^jG) • 
In summary, we have 


(117) 


/*1/G ma x r /*1/G ma x ^ ^ 

t(p)= / \s ph (e,p)\ 2 ddd(i/G) + x / V 

^0 Jo d J 0 • 1 


daj 

d(l/G) 


2 

d(l/G). 
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l/G 


da i 

OL 2 

da 2 

as 

da 3 

9(1 /G) 

9(1 /G) 

9(1 /G) 

0.00 

0.702988 

0.009776 

0.260661 

-0.017374 

0.833321 

-0.000611 

0.05 

0.705833 

-0.009915 

0.253348 

-0.046566 

0.832408 

-0.036116 

0.10 

0.704294 

-0.053006 

0.251395 

-0.029803 

0.829828 

-0.066179 

0.15 

0.700617 

-0.097783 

0.250099 

-0.016222 

0.825956 

-0.087744 

0.20 

0.694664 

-0.144215 

0.249306 

-0.010052 

0.821312 

-0.096545 

0.25 

0.686959 

-0.169986 

0.247309 

-0.061204 

0.817120 

-0.066627 

0.30 

0.677167 

-0.227359 

0.243807 

-0.072388 

0.815138 

-0.008931 

0.35 

0.664000 

-0.306018 

0.239969 

-0.074632 

0.816970 

0.085964 

0.40 

0.645668 

-0.434744 

0.237317 

-0.026502 

0.823706 

0.183724 

Table 1. Coefficients two-dimensional IOFD 


l/G 

Q(i 

da i 

0(2 

da 2 

0(3 

da 3 

0(4 

da 4 

0( 5 

da 5 

9(1/0) 

9(1 /G) 

9(1 /G) 

9(1/G) 

9(1/0) 

0.0000 

0.635413 

-0.000228 

0.210638 

0.016303 

0.172254 

-0.014072 

0.710633 

-0.006278 

0.245303 

0.019576 

0.0500 

0.635102 

-0.015578 

0.210152 

-0.023424 

0.171912 

-0.005802 

0.709821 

-0.047764 

0.245148 

0.021398 

0.1000 

0.634166 

-0.034804 

0.208167 

-0.043396 

0.171146 

-0.012462 

0.707374 

-0.070981 

0.244762 

0.007493 

0.1500 

0.632093 

-0.054496 

0.205348 

-0.065935 

0.170031 

-0.022145 

0.703359 

-0.088202 

0.245160 

0.009937 

0.2000 

0.628341 

-0.103457 

0.201605 

-0.069385 

0.169740 

0.001893 

0.698813 

-0.092327 

0.245687 

0.012201 

0.2500 

0.622526 

-0.133896 

0.197423 

-0.098212 

0.169475 

-0.002559 

0.694726 

-0.066617 

0.246454 

0.016791 

0.3000 

0.614611 

-0.183988 

0.192414 

-0.115398 

0.168690 

-0.005589 

0.692615 

-0.011177 

0.247743 

0.029213 

0.3500 

0.603680 

-0.255991 

0.186819 

-0.120930 

0.167581 

-0.015564 

0.694109 

0.077605 

0.250098 

0.059733 

0.4000 

0.588498 

-0.356326 

0.180737 

-0.132266 

0.166640 

-0.001852 

0.700902 

0.199685 

0.254352 

0.106049 


Table 2. Coefficients three-dimensional IOFD 


This integral is discretized, using a weighted sum of regularly sampled contributions. By 
n a we denote the number of angles to discretize the integration over the sphere and by 
tiq the number of choices for the parameter l/G. In two dimensions we chose = 20, in 
three dimensions = 200. The l/G axis was discretized in steps of 0.01 in the integral 
(117). The value of A = 10 -12 was used and control nodes where chosen in the interval 
[0,0.4] with distance 0.05. 

The objective functional was minimized using the Matlab function lsqnonlin, aimed 
particularly at least-squares problems. We found that the optimization problem using the 
least squares objective functional converges better than other types of objective function¬ 
als, such as a sup-norm. This was done in three steps. First the control values for l/G in 
[0,0.2] were determined, then for l/G in [0.1, 0.3] keep the values for l/G < 0.1 equal to 
those already obtained, and then for l/G in [0.2, 0.4] keeping those for l/G < 0.2 already 
obtained. In this way somewhat better result were obtained than when minimization was 
done directly for G G [0,0.4]. The parameters were determined heuristically, in such a way 
that increasing the number of discretization points would not yield substantial improve¬ 
ments. The phase slowness errors for l/G G [0,0.1] were most sensitive to details of the 
method. As the phase speed errors are very small in this parameter range we have not 
explored this further. 

The results of the optimization for the two- and three-dimensional case are given in 
Tables [l] and [2j The phase slowness errors as a function of l/G (maximum over angle) are 
given in Figure [3j together with those of the IOFD and CH06 methods. 


5.2. Amplitude correction operators. Here we will construct the difference operators 
Q. We will use the second order discretizations of the identity defined in ( |99| ). Here we 
will denote the coefficients of this family by /3y, which will be functions of = l/G. Thus 
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phase slowness errors of IOFD, QS-FEM, SUT and CH06 



Figure 3. Phase slowness errors for the IOFD method compared to QS- 
FEM and the sixth order method of m and m- 


the functions g 1 associated with Q, see Q are given by 


(118) 


in 3-D and by 


— 


Pi 

for 7 € {- 1 , 0 , l} 3 , 7 = 0 

fa 

§ 

for 7 € {- 1 , 0 , l} 3 , 7 = 1 

P3 

12 

for 7 € {- 1 , 0 , l} 3 , 7 = 2 

1-01-02-03 

8 

for 7 € {-1,0, l} 3 , | 7 | = 3, 


(119) 


— 



for 7 € {- 1 , 0 , l} 2 , | 7 | = 0 
for 7 € {- 1 , 0 , l} 2 , | 7 | = 1 
for 7 € {- 1 , 0 , l} 2 , | 7 | = 2 , 

in 2 -D, where /3j = /3j(l/G). The /3j(l/G) will be defined by Hermite interpolation from 
control values similarly as we did for the aj(l/G). The control values are chosen to 
minimize a discrete approximation of the integral 

2 


( 120 ) 


»1/G„ 


'Qd 


Q(0- 


m/omi 


d6d(l/G). 


Z=us ph (0)0 


This integral is discretized in the same way as in the previous subsection. This results in 
a linearly constrained linear least squares problem which is easy to solve in Mat lab. The 
resulting coefficients are given in Tables [3] and [4] below. The maximum over angle of the 

varied between around 10 -8 for 1/G = 0.05 and 10 -2 for 1 /G 

£=ws ph (6)6 

= 0.4. 

5.3. Comparison of phase slowness errors. The following conclusions can be drawn 
from the data in Figures [lj [2] and [3j First the QS-FEM method of |3] (in two dimensions) 
and the IOFD method developed here (in two and three dimensions) perform remarkably 
well considering their small stencils. They provides a substantial improvement, roughly a 


error 
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l/G 

pi 

a/3i 

@2 

dfo 

3(1 /G) 

3(1 /G) 

0.00 

0.872589 

-0.115476 

0.088139 

0.232493 

0.05 

0.870989 

-0.080799 

0.089351 

0.080994 

0.10 

0.866560 

-0.122182 

0.092018 

0.075452 

0.15 

0.858994 

-0.189920 

0.096178 

0.106183 

0.20 

0.847495 

-0.277477 

0.102309 

0.147420 

0.25 

0.830913 

-0.394429 

0.110797 

0.198380 

0.30 

0.807375 

-0.559277 

0.122158 

0.261263 

0.35 

0.773715 

-0.806746 

0.137030 

0.337561 

0.40 

0.724163 

-1.211119 

0.155971 

0.420753 

Table 3. Coefficients amplitude correction operator Q in 2-D 


l/G 

Pi 

d/3! 

P2 

d/3 2 

Ps 

d/3 3 

3(1 /G) 

3(1 /G) 

3(1 /G) 

0.0000 

0.806683 

0.002423 

0.193113 

-0.002685 

-0.056266 

-0.002551 

0.0500 

0.832963 

-0.081724 

0.114016 

0.032813 

0.020075 

0.058590 

0.1000 

0.841034 

-0.130484 

0.076623 

0.029868 

0.061360 

0.078398 

0.1500 

0.833587 

-0.231333 

0.076280 

0.129614 

0.067935 

0.024410 

0.2000 

0.821230 

-0.304691 

0.078943 

0.086321 

0.074389 

0.130587 

0.2500 

0.803736 

-0.416375 

0.081855 

0.072002 

0.084073 

0.220607 

0.3000 

0.779384 

-0.573760 

0.084646 

0.054207 

0.098065 

0.329810 

0.3500 

0.745468 

-0.801027 

0.086156 

0.004734 

0.118341 

0.486328 

0.4000 

0.697405 

-1.148951 

0.083351 

-0.136764 

0.148391 

0.732785 


Table 4. Coefficients amplitude correction operator Q in 3-D 


factor 20, in phase errors compared to the compact sixth order scheme SUT and CH06 of 
m and [29], which in turn are better than other alternatives. For higher order FD and 
FE methods, as can be expected, the error becomes small if both the number of points 
per wavelength and the order N become large, however this effect sets in quite late, e.g. 
at eight points per wavelength and N — 8 the relative phase slowness errors of the finite 
element method are roughly equal to those of QS-FEM and IOFD. 

Next we discuss how much accuracy might be needed, and in how far the improvements 
will make a difference in simulations. In view of © it is not unreasonable to require at 
least that 5 p h < 0.01^. In a regime of wave propagation over several hundreds wavelengths, 
using a mesh with five points per wavelength, from the methods considered only QS-FEM 
and IOFD satisfy this. At six points per wavelength the CH06 method is near this 
bound while FE8 (which is much more expensive) also qualifies. So in these situations the 
improved phase slowness accuracy obtained by using QS-FEM or IOFD can be expected 
to have some impact in terms of lower cost compared to FE8 and in terms of improved 
accuracy compared to CH06 and other compact finite difference methods. The latter will 
be confirmed in the examples in the next section. 

6 . Numerical examples 

In this section we present two numerical experiments, first in a constant medium, and 
then in a smoothly varying medium. We will present two-dimensional examples with large 
domain sizes on the order of hundreds of wavelengths. 

As mentioned, phase slowness errors typically lead to phase shift errors in the solutions. 
Considering wave propagation over 500 wavelengths as an example, it follows from © and 
the surrounding discussion that these phase shifts errors for IOFD should be negligibly 
small for meshes with five or six points per wavelength, and still quite small for four 
and three points per wavelength. For other methods these errors should show up much 
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0 0.2 0.4 0.6 

theta 


Figure 4. Phase shift errors are easily observed by plotting a 45-degree 
part of an annulus, see figure (b). Cubic spline interpolation is applied to 
map the data of figure (a) to polar coordinates. 

stronger. In our first example we will verify this numerically, assuming a constant velocity 
model. 

To simulate a point source at a given grid point, we will simply use a discrete 5-function. 
An unbounded domain is simulated by adding a damping layer around the domain of 
interest, with a nonzero imaginary contribution to k that quadratically increases from 
the boundary of the domain of interest^} The discrete system of equations is formed 
using a Matlab code written for this purpose and then either solved directly or, for the 
larger examples, exported to disk. In the latter case, the resulting linear systems are then 
solved using the MUMPS parallel direct solver [2] on a few nodes of the Lisa cluster of 
surfsara (www.surfsara.nl). This system contains 32 parallel nodes with each two intel 
Xeon processors E5-2650 v2 running at 2.60 GHz and 64 GB memory, connected by 
Mellanox FDR Infiniband. In the examples in of section between 1 and 4 nodes were used 
in parallel. 

To easily observe the absence or presence of the phase shifts, we plot the resulting 
wave field on a 45 degree segment of an annulus, with the radial coordinate varying on an 
interval of about a wavelength. The location where the real part is minimal, according to 
the exact solution, is indicated by a line that is plotted. The transformation of the field 
to polar coordinates is done by using cubic interpolation from the numerical solution on a 
Cartesian mesh. Schematically this is displayed in Figure [4j where part (b) of the figure 
is a plot in polar coordinates of the indicated region of part (a). 

The results from the computations are displayed in Figure [5j Part (a) shows that for 
second order finite differences at 10 points per wavelength (ppw) a clearly visible phase 
shift already occurs after 20 wavelengths. In (b) we see that for the JSS method a clearly 
visible phase shift occurs after 50 wavelengths. In (c), (d) and (e) we investigate the sixth 
order method CH06 of |29j at 6, 5 and 4 ppw. (We have chosen one of the higher order 


2 In a 1-D damped Helmholtz equation — — k 2 u with k constant, k = a + i/3 solutions decay as 

u — e ' l ( cx +' l P) x ' if ^ varies slowly the damping becomes proportional to e~ 1 ^00 dx m The quadratic profile 
is chosen such that e~^^ dx is on the order of 0.001 to 0.01. Reflected waves pass twice through the 
damping layer. Unfortunately reflections occur due to the medium variations. To make these small Im (k) 
must increase slowly and these layers must be quite thick. In our experiments we used on the order of 5 
to 10 wavelengths. 
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(a) 


(f) 


polar plot FD2 at 10 ppw 
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polar plot IOFD at 6 ppw 
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(b) 


polar plot JSS at 6 ppw 
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( c ) 


polar plot CH06 at 6 ppw 
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( e ) 

polar plot CH06 at 4 ppw 
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(g) 


polar plot IOFD at 5 ppw 
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(h) 
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polar plot IOFD at 4 ppw 
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(i) 


polar plot IOFD at 3 ppw 
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(j) 


polar plot IOFD at 2.5 ppw 


— 


0.2 0.4 0.6 

theta 


Figure 5. Plots of numerical solutions over a 45 degree part of an annulus 
for several numerical methods, (a) FD2 at 10 ppw; (b) JSS at 6 ppw; (c), 
(d), (e) sixth order method of [29] at 6, 5 and 4 ppw; (f)-(j) IOFD method 
at 6, 5, 4, 3 and 2.5 ppw. 


methods). At 6, 5 and ppw the maximum phase errors at 500 wavelengths are 0.27, 0.84 
and 7r radians respectively and the associated phase shifts are increasingly visible in the 
pictures. In parts (e) to (i) we plot results for the IOFD method at 6, 5, 4, 3 and 2.5 
ppw. At 6, 5 and 4 ppw the maximum phase shifts are 0.0065, 0.020 and 0.089 radians 
respectively, i.e. considerably smaller than observed for CH06. At 3 ppw the phase shift 
after 500 wavelengths is clearly visible, only at 2.5 points per wavelength does it become 
large and in this case the field is plotted at 100 instead of less than 500 wavelengths from 
the source. 

In Figure [6] we plot the amplitude errors for IOFD at 3 and 4 ppw. For more than 4 
ppw they were increasingly small. 

In our second example k is variable. To avoid that errors due to the discretization 
of the velocity model become dominant we use use a smoothly varying velocity model, 
namely a smoothed Marmousi model. In this example we will compare a solution with 
IOFD using a minimum of six points per wavelength with a fourth order finite element 
solution using twice as many grid points in each direction. In these examples the right 
hand side was a point source and the linear systems were again solved with MUMPS. In 
case of variable coefficients we assumed that k(x) is defined on the cell centers. The values 
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(a) 

amplitude errors IOFD at 4 ppw xio 4 



0 0.2 0.4 0.6 



0 0.2 0.4 0.6 


Figure 6. Amplitude errors relative to exact solution for IOFD at (a) 4 
ppw (b) 3 ppw. 


smoothed Marmousi velocity model 



0 1000 2000 3000 4000 5000 6000 7000 8000 9000 

x(m) 


Figure 7. Smoothed Marmousi velocity model 


of f 1 {kh(x )) (see (54)) at other points were obtained using linear interpolation from the 


values of f 1 {hk{x )) at the cell centers. 

The velocity model is given in Figure [7| It is obtained from the Marmousi model by 
convolving along both of the axes with a cos square pulse of width 160 meter. We will 
give results for 50 and 100 Hz. A solution for the first case is given in Figure [8j Figure [9] 
contains four plots. The top plots are reference amplitudes for obtaining relative errors 
and the bottom two plots are relative errors with respect to the reference values. In both 
cases we give results for 50 and for 100 Hz. The reference value is a local average of the 
absolute value of the solution over a square of about 2 by 2 wavelengths. This is done 
because the solutions themselves contain nodal points from interfering waves, where the 
amplitude is very small, and are hence not directly suitable as reference value. Very small 
relative errors are obtained (except directly at the source point), ranging from less than 
0.01 over most of the domain to 0.05 or 0.08 at isolated spots where the absolute amplitude 
is small. 


7. Application in multigrid based solvers 

The last few years there have been several interesting developments in multigrid methods 
for Helmholtz equations. Different two-grid methods with inexact coarse level solvers 
have been studied in [5] and [25j. In [5] a number of iterations of shifted Laplacian 
preconditioned Krylov solver nn is used as coarse level solver. The method of is based 
on the multigrid method in |[26j with a double sweep domain decomposition preconditioner 
m as coarse level solver. The multigrid method with exact coarse level solver was studied 
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solution at 50 Hz 



0 1000 2000 3000 4000 5000 6000 7000 8000 9000 

x(m) 

Figure 8. Solution from a point source at 50 Hz 



Figure 9. Reference values for 50 Hz (a) and for 100 Hz (b). Relative 
errors in the solutions for 50 Hz (c) and for 100 Hz (d). 


in [26] . There it was shown that the convergence can be strongly improved when phase 
slowness differences between the fine and coarse scale operators are minimized. For this 
purpose, optimized finite differences were used at the coarse level, and good convergence 
was obtained for meshes with downto three points per wavelength at the coarse level. For 
standard choices of the coarse level discretization it was found that about 10 points per 
wavelength at the coarse level were needed to have good convergence. 

In [26] standard second order finite differences were used as the fine level. Because 
of the relatively large phase slowness errors of this method, the coarse level optimized 
finite difference method had to be constructed specifically to match the phase slownesss 
of second order finite differences, instead of matching the true phase slowness. A better 
choice is to use method with small phase slowness errors at the fine level and at the coarse 
levels. Here we will use IOFD at all levels. These experiments do not involve the operator 
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Figure 10. Velocity models for the 2-D two-grid experiments, (a) Mar- 
mousi (b) 2-D slice of the SEG-EAGE salt model. 


Q. The operator P is used directly as coarse level discretization and at the fine level we 


are only interested in solving equation (111). 


In the first set of computational results of this section we will show that this results 
in good convergence of the multigrid method with exact coarse level solver. In a second 
example we will study the multigrid method with inexact coarse level solver of [25] . 

In our study of the convergence when using the exact coarse level solver we are again 
interested in examples with wave propagation over hundreds of wavelengths. Therefore, 
these experiments are done in two dimensions. For background on multigrid methods, 
see [28] . As in [26J most of the components of the multigrid method are standard. Full 
weighting restriction and prolongation operators are used. As smoother, an cj-Jacobi 
method is used. We found that co — 0.7 and v — 4 (the number of pre- and postsmoothing 
parameters) are good choices of parameters. In these experiments we used a conventional 
absorbing boundary layer to simulate an unbounded domain. 

We studied the convergence as a function of the number of points per wavelength for 
three velocity models: A constant model, the Marmousi model and a slice of the 3-D SEG- 
EAGE salt model. The latter two models are displayed in Figure [lO] The parameters of 


the examples and the observed number of iterations to reduce the residual by 10 b are given 
in [5| At downto three points per wavelength the method behaved well. At 2.5 ppw coarse 
level the method still converged, but the number of iterations increased substantially, and 
also became more sensitive to the problem size (which was apparent from smaller scale 
experiments not included in the table). 

Note that the application in multigrid methods is quite different from the application 
as fine level discretization. The method is used at coarser meshes (at three points per 
wavelength the direct application will in general lead to too large errors). Also, multi¬ 
grid solvers using IOFD at coarse levels may be developed for other types of fine level 
discretizations, as long as they use a regular mesh. 

We now turn to a multigrid method with an inexact coarse level solver. Such methods 
are used because in three dimensions it is often too expensive to compute the exact solu¬ 
tion. These methods are currently some of the fastest solvers for large problems that are 
in the literature [a E5]. 

Because we are interested in coarse meshes, such as six points per wavelength based on 
the previous examples, it is a priori not clear that the above mentioned solvers perform 
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constant 

2400 x 2400 

Marmousi 
4600 x 750 

salt model 
2700 x 836 

ppw 

freq 

its 

freq 

its 

freq 

its 

5 

480 

29 

150 

23 

60 

18 

6 

400 

8 

125 

11 

50 

8 

7 

342.9 

6 

107.1 

9 

42.9 

7 

8 

300 

5 

93.8 

8 

37.5 

6 

9 

266.7 

5 

83.3 

7 

33.3 

6 

10 

240 

4 

75 

6 

30 

5 


Table 5. Iterations required for a two-grid method using IOFD discretiza¬ 
tion at the fine and coarse level as a function of the number of points per 
wavelength (ppw) 


well. Like many solvers in the literature, they were tested for problems with at least 
ten mesh points per wavelength. They cannot be assumed to converge as well for larger 
frequencies, because multigrid convergence depends on frequency, and the same is true for 
the shifted Laplacian preconditioner [9J. For the double sweep domain decomposition it is 
unclear how the frequency affects the convergence, but a priori it also cannot be assumed 
to be independent of the frequency. 

This raises the question whether we can actually obtain a gain in efficiency by going 
to coarser meshes. The purpose of the next example is to show that this indeed the case, 
and to generally show that IOFD can perform well with the solver of [25] . 

In the following example we will test the method of [22], which is a two-grid method 
using an inexact coarse level solver given by a double sweep domain decomposition pre¬ 
conditioner (see [211). The method is modified to use IOFD at both the fine and the coarse 
levels of the two-grid method. We will take the SEG-EAGE Salt Model as an example, 
similarly as in [25] . In addition to changing the discretization method we will increase the 
frequency by a factor |, so that a minimum of six points per wavelength is used, a regime 
which has not been tested before for this method. If convergence and cost per degree of 
freedom would stay constant, there would be an improvement in the cost by a factor of 
over (g) r^> 4.62 (more than this because cost grows somewhat faster than linear with 

problem size). 

The original SEG-EAGE salt model is of size 13500 x 13500 x 4200 meter, discretized 
with 20 m grid spacing. We apply the method just described to solve the Helmholtz 
equation with this velocity model and random or point sources as right hand sides at four 
different frequencies from 6.25 to 12.5 Hz. Slices of the model are displayed in Figure [Tl] 
Parameters in the two-grid method are v — 3 for the number of pre- and postsmoothing 
steps and uo = 0.65 in the c^-Jacobi method. Computations were done on the Lisa cluster 
at Surfsara, described already in section [6j using the implementation described in [25] . A 
maximum of 16 nodes were used in parallel for these computations. 

The algorithm is set up to solve for multiple right hand sides simultaneously. In the 
table of results, the computation time per right hand side is given. In Table [6] some 
parameters are given, together with the computation time and iteration count to reduce 
the residual by 10 -6 . As illustration, plots of a solution are given in Figure [l2| It can be 
observed that the cost increases very little compared to the results of [[25] . even though 
frequencies are increased by a factor 5/3. Some increase in cost can be expected, because 
the discrete Helmholtz operator using second order finite differences is cheaper to apply 
than the one using a compact 27-point stencil. Hence reducing the number of points per 
wavelength in the mesh can indeed lead to corresponding savings in computation time. 
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Figure 11. SEG-EAGE salt velocity model: (a) (x,z) slice at y — 6740 
m (b) (y, z) slice at x — 6740 m. 


(a) 


(b) 
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Figure 12. Solution to the Helmholtz equation at 12.5 Hz: (a) (x,z) slice 
at y — 6740 m (b) (y, z) slice at x — 6740 m. 


frequency 

6.25 

7.87 

9.91 

12.5 

size 

338x338x106 

426x426x132 

536x536x166 

676x676x210 

#dof 

1.3 • 10 7 

2.5 • 10 7 

5.0 • 10 7 

1.0- 10 8 

cores 

32 

64 

128 

256 

# of rhs. 

1 

2 

4 

8 

iterations 

12 

12 

13 

15 

computation time/rhs(s) 

26 

35 

45 

73 


Table 6. Computation times and iteration counts for the SEG-EAGE Salt 
Model example. 


As mentioned, the methods of [5J and [25j are some of the fastest currently in the 
literature. Comparing with these results we see a significant improvement. For example 
in [5] the SEG-EAGE salt model problem was solved at 10 Hz in 270 seconds on 256 cores 
of an IBM BG/P machine (with the residual reduced by a factor 10 5 instead of 10 6 in our 
case). Here we solve the problem at 9.91 Hz using 128 cores in 45 seconds per right hand 
sides (179 seconds for four right hand sides), a clear improvement]^} 

8. Discussion 

Here we summarize some of the conclusions and further discuss the results. 

Using the results presented one can make a case for the use of coarse meshes using a 
minimum of five or six points per wavelength in time harmonic wave simulations in case 
k is smooth. This idea is not new, in the exploration geophysics community it appears to 
be quite common. However, we found that the methods that have been proposed for this 
purpose in m and m can be expected to give substantial phase errors in simulations 
of large distance wave propagation. By using the new IOFD method (in two or three 

^On the other hand the method of [5] uses less memory and has been applied to larger examples than 
we have shown here. Furthermore no full comparison including accuracy was made. 
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dimensions) or the QS-FEM method (in two dimensions only), phase errors can be made 
much smaller. 

In k(x) has strong gradients, one can expect that, at least locally where \7k is large, finer 
meshes and/or different discretizations are needed to obtain accurate solutions. Large gra¬ 
dients lead to reflections. Typically finer meshes are needed to model these accurately. For 
one reason this is because finer discretizations of k are needed, since linearized scattering 
theory shows that reflected waves are associated with perturbations in the medium veloc¬ 
ity with wave vectors of length up to 2k (where here k refers to the background velocity 
around which the linearization is applied). However, in this case the multigrid approach 
discussed in section [7] can still be useful. It has been applied in successfully in examples 
with discontinuities. This suggests to do further research on multigrid approaches with 
compact finite difference method at the coarse level and other discretizations at the fine 
level, including methods with local refinement. A similar argument can be held for the 
discretization of the right hand side / in the equation (—A — k(x) 2 )u = /. For rapidly 
varying functions /, finer meshes may be needed at least locally where the rapid variations 
occur. 

When applied in inversion algorithms IOFD and QS-FEM are somewhat more com¬ 
plicated than the methods of m and m, because the operatore depends in a more 
complicated fashion on the coefficients, which means it is more complicated to compute 
the derivative of the finite difference operator with respect to the medium coefficients. 
Due to the use of Hermitian interpolation, these derivatives are however continuous for 
our IOFD method. 
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Appendix A. Phase slowness computations for finite element methods 


In a periodically repeating setting, which is the case for finite elements with N > 2, 
the situation is somewhat more complicated. The symmetry property (109) only holds for 
p, g, r divisible by TV. For such operators we consider the Bloch waves 

7/7 — pif ,' x l , m , n lh 

^hm^n c u l,m,n 


( 121 ) 


where vi, m ,n is periodic with shifts (pTV, qN , rTV), p, g, r integers. For given £, the action of 
an operator A with these symmetries is given by an TV 3 x TV 3 matrix acting on the vi, m , n 
(in three dimensions) for 0 < Z, m, n < TV — 1. We need to find the vectors £ for which 
there is a zero eigenvalue. However not all zero eigenvectors correspond to plane waves 
with wave vectors £, because of the presence of v^ mn . In general v^ m ^ n can correspond to 
a linear combination of plane waves with wave vectors given by (p-j^, r Wh)^ where 
p, g, r are integers. Assuming that some £ corresponds to a simple zero eigenvalue and 
that ui^ n is close to a plane wave (which is often the case because the eigenfunctions 
of the continuous operator are plane waves and the operator A is a good approximation 
of the continuous operator), these integers p, g, r can be determined modulo TV, and the 
wave vector associated with an element of the zero set of A can be determined. 

In the computations we will take a somewhat different approach. We will compute 
all eigenvalues, and then only consider the one whose eigenvector vi^ n is most closely 
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correlated with (has the in absolute value largest inner product with) the constant function 
— 1- We will say we have found a phase velocity vector at some £ if this eigenvalue 
is zero. This approach has some limitations, but a more extensive study of this topic falls 
outside the scope of this paper. For standard finite elements and h such that the mesh 
has more than four points per wavelength this appeared to be sufficient. 



